Predicting the potential distribution of Dactylorhiza hatagirea (D. Don) Soo-an important medicinal orchid in the West Himalaya, under multiple climate change scenarios

Climate variability coupled with anthropogenic pressures is the most critical driver in the Himalayan region for forest ecosystem vulnerability. Dactylorhiza hatagirea (D.Don) Soo is an important yet highly threatened medicinal orchid from the Himalayan region. Poor regenerative power and growing demand have resulted in the steep decline of its natural habitats populations. The present study aims to identify the habitat suitability of D. hatagirea in the Western Himalaya using the maximum entropy model (MaxEnt). The community climate system model (CCSM ver. 4) based on representative concentration pathways (RCPs) was used to determine suitable future areas. Sixteen least correlated (< 0.8) bioclimatic, topographical and geomorphic variables were used to construct the species climatic niche. The dominant contributing variables were elevation (34.85%) followed by precipitation of the coldest quarter (23.04%), soil type (8.77%), land use land cover (8.26%), mean annual temperature (5.51%), and temperature seasonality (5.11%). Compared to the present distribution, habitat suitability under future projection, i.e., RCP 4.5 and RCP 8.5 (2050 and 2070), was found to shift to higher elevation towards the northwest direction, while lower altitudes will invariably be less suitable. Further, as compared to the current distribution, the climatic niche space of the species is expected to expand in between11.41–22.13% in the near future. High habitats suitability areas are mainly concentrated in the forest range like Dharchula and Munsyari range, Pindar valley, Kedarnath Wildlife Sanctuary, West of Nanda Devi Biosphere Reserve, and Uttarkashi forest division. The present study delineated the fundamental niche baseline map of D. hatagirea in the Western Himalayas and highlighted regions/areas where conservation and management strategies should be intensified in the next 50 years. In addition, as the species is commercially exploited illegally, the information gathered is essential for conservationists and planners who protect the species at the regional levels.


Introduction
Climate is an important ecological and abiotic factor affecting species potential geographic distribution and ecological niche space [1]. The minimal changes in species bioclimatic envelope are thought to have considerable impacts on the plant-pollinator relationship, seed set, and regeneration status. It is expected that species may no longer adapt to a set of environmental conditions to facilitate further expansion [2]. Therefore, the species must either cope with the prevailing ecological conditions or colonize to sustain or become extinct [3]. This has led to a growing interest in developing and scaling up prioritization strategies for such species to ensure the highest conservation gains [4].
The Himalayan regions are an assemblage of biodiversity hotspots [5]. However, the ongoing disturbance exacerbated by climate change, habitat fragmentation, invasions by alien species, grazing and trampling, overexploitation, and excessive consumption of natural resources has altered the structural and functional integrity of the various Himalayan ecosystems [6]. Besides these, climate variability, land-use change, and rural migration are key contributors to biodiversity loss in the region [7]. In the last few decades, the region saw cascading effects of climate variability mainly due to increased greenhouse gases concentration. It is believed that the rate of global warming in the Himalayas is much higher than the global average. For instance, in the last 100 years, the global average temperature rise was 0.74˚C [8]. However, in the Himalayas, a 1.5˚C temperature increase was documented in the final quarter (i.e.,  of the twentieth century [9], with warming potentially reaching 5˚C by the end of the twenty-first century [9,10]. This rise is alarming because Himalayan floras are alienated to specific elevation gradient/microhabitat conditions [3,11]. The shift in the climatic envelope is expected to bring significant change in the resident species habitat conditions, leading to changes in species richness, population structure, and those unable to cope are likely to face local extinction [1][2][3]. Furthermore, the recent upsurge in herbal or its derived products across the globe has led to uncontrolled abusive practice; thus, the natural stock of these plants is under tremendous pressure. In the case of the Indian Himalayan Region (IHR), a considerable number (1748 species) of medicinal and aromatic plants (MAPs) are reported, with 31% of them being native, whereas 15.5% are endemic and threatened [12]. The high potential instability and inherent vulnerability make the region one of the most ecologically fragile bio-geographic zones [13]. Other challenges on the MAPs include low population size, habitat specificity, genetic bottleneck effect, narrow distributional ranges, and heavy livestock grazing [14]. The literature on these threatened plants is fragmentary or limited to specific geographic pockets [15]. In the above context, it is obligatory to make a conservation framework encircling species habitat restoration and promote cultivation, thus, reducing pressure on the wild populations.
The development of statistical modelling and geospatial technology in predicting suitable habitat distribution has gained immense popularity. However, such information is at an initial phase for the Himalayan MAPs [6]. The use of geospatial technology could add an advantage as obtaining specific distribution maps for such species is difficult and often requires intensive surveys [16]. The difficulty level becomes amplified in the Himalayan region where the working conditions are not conducive for the survey, i.e., inaccessible and difficult topography perplexed with hostile conditions. Therefore, estimating current plants distribution and identifying important climatic refugia will help predict future distribution patterns and reveal regions with high extinction rates. At present, the common method to study potentially species distribution and environmentally suitable habitats is to use species distribution models (SDMs) [17]. SDM has made it possible to analyze the environmental drivers of species distributions and project a species realized niche into a geographic area [18]. Of many SDM algorithm methods, MaxEnt has proved decisive when modeling rare species with narrow ranges [19]. MaxEnt modelling is a robust computational algorithm that works on the backdrop of species presence points and rasterized environmental data. The probable ecological niche can be reconstructed using species presence data points and environmental variables/ predictors [20]. Such model-based sampling would become an important benchmark for endemics and threatened species and is a well-recognized cost-efficient method [21].
In the present study, an effort has been made to model the potential habitat distribution and effect of future climate change on D. hatagirea, a critically endangered [22,23] and endemic species of the Himalayan region. The species is a tall, slender, ground-dwelling, perennial herb with palmately lobed tuberoids that prefer to grow in a moist, mild, acidic soil environment (Fig 1A and 1B). The species has been reported from India, Afghanistan, Pakistan, Nepal, Tibet, and Bhutan [6]. In the IHR, the species is reported from Jammu and Kashmir, Ladakh, Himachal Pradesh, Uttarakhand, Sikkim, and Arunachal Pradesh at an altitudinal range of 2500 to 4500 m above sea level (asl). The estimated annual trade of the species is around 10-50 metric tons [24], with an economic value of US $ 68.88-89.54 (1US $ = Rs. 77.39) kg -1 of the dried tuber (Fig 1C). At present, the tuber of the species is destructively harvested and illegally traded; thus, it puts a stake on its future existence. Moreover, the species require specific microhabitat conditions for growth and perturbation, thus limiting the species widespread distribution. Therefore, to minimize the pressure on the wild populations, efforts are ongoing to develop and upscale the existing multiplication strategies for mass multiplication. Meanwhile, mapping and conserving the critical habitat is expected to offer a possible solution to species conservation and management. The study attempts to address the following scientific questions: (i) What is the present potential geographical distribution range of Dactylorhiza hatagirea in the Western Himalaya, India? (ii) What will be the impacts of climate variability on the future distribution of D. hatagirea using four Representative Concentration Pathways (RCPs)?, and (iii) Where are the high potential distributional areas of D. hatagirea that could be protected or could be suggested for cultivation, reintroduction/ recovery plans?. Answering these questions will help identify suitable habitats for the conservation of the species, which may help policy planers while developing strategies for its conservation.

Study area and ecological significance
The study was undertaken in Uttarakhand state (28˚43 0 to 31˚28 0 N Latitude and 77˚34 0 to 80˚03 0 E Longitude) of IHR (Fig 2). The state has a total recorded forest area (RFA) of 38,000 km 2 (71.05% of its total geographical area 53,485 km 2 ), out of which 26,547 km 2 is reserved forest, 9,885 km 2 is protected forest, and 1,568 km 2 is unclassed forests [25]. The state experience varied climates from warm dry to warm wet and with a latent cool, dry period. The state's temperature ranges from sub-zero to 43˚C [26], and average annual rainfall varies from 1093.8 mm to 1385.5 mm [27].
The present study area harbors alpine vegetation, which covers 8,524 km 2 . Of these, 4,376 km 2 is surmounted by permanent snow cover (i.e., corresponding to ca.24.11% statealpine geographical area) [28]. The alpine regions are well known for their high-value MAPs, including D. hatagirea. The region is experiencing major environmental transformation repercussions, and anthropogenic activities outnumber the natural eventualities, thus enforcing species to various threat categories [15].

Species point data
To predict species distribution, it is a prerequisite to have species presence points and environmental variables [29]. The data search was primarily made from online portals such as Global Biodiversity Information Facility [30], published literature, and herbarium consultation Botanical Survey of India, Dehradun (BSID). Data on the species were very limited, whereas herbarium records were not geo-referenced. Considering these limitations, an extensive field survey was conducted during 2016-19, and presence points were recorded. A total of 30 occurrences of the species were recorded during field surveys. A portable multi-channel Global Positioning System (Garmin) receiver with 10-20 m positional accuracy was used to record the species occurrence geo-coordinates. The coordinates were then converted to decimal degrees and used to model the species potential habitats distribution in its native range.

Data source
The climate data were downloaded from World Climate Database [31]. WorldClim provides current (baseline) and projected climate data for 2070 with a spatial resolution of 30 seconds (ca. 1 km) in GeoTIFF format. These climatic data are derivatives of maximum, minimum, and average values of monthly, quarterly, and annual temperatures and precipitation of the last 30 years, i.e., 1970-2000. Likewise, environmental variables such as soil type, soil moisture, and soil pH were downloaded from the International Soil Reference Information Centre [32], while land use land cover (LULC) from http://www.esa-landcover-cci.org/ [33] ( Table 1). Besides these, non-climatic variables, i.e., altitude, aspect derived from NASA Shuttle Radar Topographic Mission (SRTM, version 4.1) [34]. The reason behind using both the climatic and non-climatic variables is to enhance the model's predictive power as suggested for endemic plants [5,35]. Further, for future prediction studies, we used Community Climate System Model (CCSM) ver. 4 (CCSM4) that is based on the Fifth Assessment Report (AR 5) of the Intergovernmental Panel on Climate Change [36] and two contrasting Representative Concentrations Pathways (RCP 4.5 and RCP8.5) for the years 2050 and 2070. Furthermore, we assumed that the edaphic properties are expected to remain stable in the next several decades, as soil properties should not change synchronously with sudden climate change; hence the same raster layer was used in future projections.

Environmental layers and variable selection.
The model's output can be accurate, biologically meaningful, and generalized if built with predictor variables that directly impact species distribution. Strong collinearity between the variables in SDMs may cause model overfitting due to the high level of correlation among variables [37]. To avoid multi-collinearity among the 19 bioclimatic variables, highly correlated variables (r � 0.80 Pearson correlation coefficient) were eliminated from further models using ENM Tools. This reduction of predictor variables resulted in the inclusion of nine bioclimatic variables and seven environmental variables for the prediction process (S1 Table). Further, using ArcGIS 10.3, all predictor variables layers were rasterized into the same bounds, cell sizes, and coordinate system as the layer of occurrence localities. Finally, these layers were converted to ASCII format for further processing in MaxEnt.

Model parameterization. MaxEnt algorithm (MaxEnt ver. 3.4.1) [38]
for habitat distribution modelling was employed [29]. MaxEnt algorithm was chosen over other available machine learning tools owing to (i) presence of only data points of the species, (ii) works even relatively with a small number of occurrence locations and high predictive performance, (iii) can handle continuous and categorical environmental data simultaneously, (iv) analyze results in terms of percent contribution of environmental data through model output, (v) examine variables weight through jackknife method, and (vi) calibrate the model, run numerous replicates along with cross-validation, and bootstrapping to test model robustness [18,39,40]. We used 75% of the dataset for training and 25% dataset model testing in this study. For generating model robustness, the number of iterations was set to 5000, with 30 replicated model runs. The maximum background points10000 and ten percentile training presence with logistic threshold rule were applied, whereas other parameters were set to default.

Model performance and potential niche change.
To calibrate the model and validate its robustness, threshold independent receiver-operating characteristic analysis (ROC) and area under the receiver-operating characteristic curve (AUC) were tested for model precision. The AUC value varies between 0 to1. The values close to +1 indicate conformity between observations and prediction, whereas zero or less values indicate a performance no better than random [41]. Statistically, AUC values near 1 indicate very good model performance, whereas AUC values close to 0 signify complete inaccurate prediction. Model performance based on AUC values are categorized as, very good (0.95 < AUC < 1.0), good (0.9 < AUC < 0.95), fair (0.8 < AUC < 0.9), and poor (AUC < 0.8) [42,43]. In the past, several studies have suggested that the AUC values mislead the performance of predictive distribution models and reflect relative model performance [44]. Therefore, to assess the predictive success of models, sensitivity, specificity, overall accuracy, and True Skill Statistics (TSS) were calculated by a confusion matrix. Threshold-dependent TSS is considered an additional accuracy measurement that is

PLOS ONE
Potential distribution of Dactylorhiza hatagirea (D. Don) Soo not affected by prevalence as it does for the kappa coefficient and the size of the validation set [45]. It deals with sensitivity and specificity, values ranging from − 1 to + 1, where + 1 indicates perfect agreement, scores ranging from 0.6 to 0.9 specify fair to good model performance, and 0 represents a random fit [45]. For this, the output of the logistic layer derived from MaxEnt results was reclassified into a binary prediction map (unsuitable and suitable) with a threshold of 10 percentile training presence. All geographical plotting and suitable range-size estimation were conducted in ArcGIS software (version 10.3).
To identify the potential area of distribution, the distributional indices based on threshold interval classification (TIC) were categorized as highly suitable (TIC > 0.75), moderate suitable (0.50 < TIC < 0.75), least suitable (0.25 < TIC < 0.50), and unsuitable areas (TIC < 0.25). Changes in the potential niche of D. hatagirea between the current and future climatic scenarios were computed by converting ASCII output projections into raster format using ArcGIS 10.3. Simultaneously the number of cells (pixels) among the projected climatic extent was calculated using zonal statistics in spatial analyst tools in ArcGIS 10.3. The differences in the mean number of cells among four classes of potential habitats were converted to surface area (km 2 ). Finally, MaxEnt predictive maps for the current and future scenarios were related to elevation classes. This would help map habitats and contribute to species-specific interventions/ reintroduction programs.

Preliminary screening of model inputs variables
The credibility of any prediction model is dependent on input variables for species distribution modelling. Given this, sixteen predictor variables out of twenty-six variables; with correlation coefficients of � 0.8 were retained after preliminary screening and selected for further modelling (Table 1).

Model performance and variable contributions
The results obtained by an ecological model are judged for their performance based on complex algorithm tests and model validation. The threshold-independent ROC showed that the average AUC yielded satisfactory results of 0.96 (Fig 3), which falls under 'very good' (0.95 < AUC < 1.0) model performance based on Thuiller et al. (2005) [42] classification. The confusion matrices for the current prediction model calculated the model's sensitivity and specificity to be 0.79 and 0.95, respectively. With these matrices in place, the model performance (i.e., TSS) was calculated (sensitivity + specificity-1). The TSS for the current model was computed to be 0.74, which indicates that the model's overall performance was good, based on Allouche et al. (2006) [45] criteria.
The variable contributions analysis highlights; elevation had the most (34.85%) influential effect followed by precipitation of coldest quarter (Bio 19; 23.04%), soil type (8.77%), LULC (8.26%), annual mean temperature (Bio 1; 5.51%), temperature seasonality (Bio 4; 5.11%), precipitation of warmest quarter (Bio 18; 4.38%) and Slope (3.06%) ( Table 1). The variables mentioned above cumulative contributions stood at~93% to the modeled potential climatic niche of D. hatagirea. Similarly, Jackknife analysis indicates annual mean temperature (Bio 1), elevation, mean temperature of the wettest quarter (Bio 8), and precipitation of coldest quarter (Bio 19) as most important predictor variables (Fig 4). These variables provide useful and distinctive information defining the D. hatagirea distributions when used in isolation. Variables like soil type, LULC, temperature seasonality (Bio 4), and annual temperature range (Bio 7) showed considerable change and showed moderate gain when used separately (Fig 4). Furthermore, the quantitative relationship between the logistic probability and input variables are depicted as response curves (Fig 5). A geomorphic variable, such as elevation, was one of the key variables that describe the present and future spatial distributions of D. hatagirea. Response curves analysis reveals average altitude ranged from 2800 m to 4500 m, precipitation of coldest quarter (Bio 19) ranged from (150 mm to 380 mm), annual mean temperature (Bio 1) ranged in between (0 -25˚C). Likewise, precipitation of the warmest quarter (Bio 18) ranged in between (200 mm-1100 mm), soil pH (5-6), and slope angle ranged from (5˚-45˚). Thus, all the identified variables estimate the important climatic attributes that potentially influence the distribution of D. hatagirea in northwestern Himalaya, India.

Current predicted potential distribution of climatically suitable areas
Habitat suitability of D. hatagirea was determined based on threshold interval classification (TIC). Of these, maximum of 96.18% (51637 sq. km) of the geographic region was predicted to   (Fig 6). Likewise, moderate habitat suitability was located in GovindVihar National Park, Uttarkashi forest division, Kedarnath Wildlife Sanctuary, Nanda Devi Biosphere Reserve, and National Park, Pindar valley and in the forest range of Dharchula and Munsyari range.

Shifts in habitat suitability under the climate change scenarios
The final output maps (current and future scenarios) were employed to find out habitats that will remain stable, gains in habitat area, habitat loss, and unsuitable habitat as part of computing change analysis (in sq. km) (  Table 3). The contraction in habitats was mainly recorded from the low habitat suitability class. A detailed shift in degrees in a different class is tabulated in Table 4.

Discussion
The susceptibility to climate change is now reflected in spatial distribution and forest ecosystem vulnerability across the globe [46]. Climate change is projected to be a 'dominant stressor' under different climate projection models in the latter half of the 21 st century [47]. As there is no denial of the reality of climate change, attention is now being actively paid to formulate some mitigation measures to maintain the stability of an ecosystem. Species distribution models in this direction have contributed immensely by providing reliable information about the potentially suitable habitats of sensitive/ vulnerable species or communities that need priority attention. The present study investigates the habitat suitability of D. hatagirea under a changing climate scenario using a species distribution model. The study provides a detailed map of the current and future species distribution in light of climate changes (Figs 6 and 7). Dactylorhiza hatagirea distribution was mostly explained by topographical variables rather than bioclimatic variables (Table 1). Among the topographical factors, elevation (34.85%) was the most dominant contributing factor, followed by soil type (8.77%) and LULC (8.26%), while among the bioclimatic variables, precipitation of coldest quarter (Bio 19; 23.04%), annual mean temperature (Bio 1; 5.51%), and temperature seasonality (Bio 4; 5.11%) were most prevalent. These physiographic factors, along with topographical features (i.e., elevation, land use characteristics, slope angle, and aspect), and bioclimatic parameters, are reported to have pronounced effects on the pattern of species distribution and community structure in the alpine meadows [35,48]. Such factors play a key role in the alpine biodiversity where the species are skewed more towards particular habitats (i.e., moist and marshy habitat) than open grasslands or rugged terrains [14,49]. Previous findings are in line with this study, wherein altitude and bioclimatic variables (temperature and precipitation) were reported to play a major role in the distribution and population structure of D. hatagirea [50][51][52][53]. Among these, Thakur et al. (2021) [52] reported     [3], using ensemble species distribution modelling (eSDM), reported elevation (30.97%) as the major dominant contributing variable, while amongst the bioclimatic variables, the mean temperature of the wettest quarter (Bio 8; 24.69%) and annual precipitation (Bio12; 21.11%) were the key contributing variables. Given the understanding from the above and several other studies on terrestrial orchids, other than biotic elements, temperature (of the growing season) and precipitation strongly modulate or affect their distribution. Thus, any significant climate change can be envisaged to have a magnified impact on these species overall distribution and growth performance. An observational study in the laboratory conditions noted maintaining the live germplasm of the species at controlled temperature chambers, i.e., 20˚C, 30˚C and in glasshouse condition (>35˚C), revealed better growth performance and flower development at 20˚C only. In comparison, growth ceased after the initial leaf development stage and later perished in temperature above 35˚C, suggesting species narrow temperature tolerance regime for growth and perpetuation. Similar findings are reported elsewhere, where an excessive rise in temperatures was reported to affect the vegetative growth flower bud differentiation negatively and preclude regeneration of Paeonia delavayi [57], Camellia sinensis [37], Dalbergia cultrata [58], Rosa arabica [59], Hippophae salicifolia [60], Fritillaria cirrhosa [61], Rhododendron niveum [62], and Ilex khasiana [43]. In addition to the above, soil type, moisture ratio, and soil pH [49] are the other vital decisive variables that posture or play a role in shaping the distribution of D. hatagirea (Table 3). In a study by Thakur et al. (2021) [52] revealed that the populations of D. hatagirea flourished in soils rich in organic matter (i.e., loamy sand to sandy loam) and had adequate soil moisture  [49]. Much of these habitats had slope steepness angles ranging from gentle (14˚) to moderate slope (41˚) angle, while north/ north-east aspects were most prevalent [49]. These variables were identified as factors of importance in our study, considering their role in posturing the vegetation reserve in the alpine habitats. Therefore, their role cannot be ignored in the Himalayas context, where the influence of the microclimatic conditions edaphic factors largely prevail [63].
In recent years, rapid climate changes in the Himalayas have resulted in distributional changes for a wide range of taxa [64]. The aftermath of these events has seriously impacted the geographical distribution of species, with reports of some species migrating to higher elevations [1,57,[65][66][67][68]. The alpine species in the Himalayas are habitat-specific and have a narrow distributional range; therefore, they are more vulnerable to extinction. The disproportionate effects on these species migrating to high latitude or elevation are attributed to climate change [65]. Our prediction showed that the shift in suitable habitat distribution to high elevations would gradually become more significant using the climate projection model. Model predicted climatically suitable habitat would expand under the RCP 4.5 and RCP 8.5 climate scenario toward north/ northeast direction and would become invariably less suitable at the lower altitudes (Figs  6 and 7). Further, we speculate that the north/ northeastward shift observed in our study could also be attributed to; North-facing aspects in the mountains are associated with higher soil nutrient content, higher biomass, coverage, height, species diversity than south-facing slopes [49,69,70]. Climate warming and shifting of species towards north/ northeast in the Himalayas is reported [3,35,71,72]. In Nepal, Shrestha et al. (2021) [53] projected that the target species would be elevated to 5000 m in the future, a substantial change when compared with the present distribution at 4000 m. A similar prediction for the species has been made elsewhere [3,50].
Likewise, the change analysis, i.e., stable area, habitat loss, and habitat gain area, was carried out for RCP 4.5 and RCP 8.5 (Table 3). The area under stable habitat under RCP 4.5 (2050) was invariably higher RCP 8.5 (2050 and 2070); suggesting changes in climatic parameters and land use characteristics over higher emission rate would negatively impact the suitability of the habitat. Other than these, model prediction advocates climatic changes will also bring new areas (habitat gain) under habitat suitability class, with the highest gain observed under lower emission rate, i.e., RCP 4.5 than RCP 8.5. The study prospects the RCP4.5 scenario to be more favorable for D. hatagirea than RCP 8.5 in the northwestern part of IHR, thus providing an expansion scope. Meanwhile, warming of lower elevation or areas where the species is currently found explains the climatic niche loss in the future. The result obtained by overlaying the projected layer in the current and future climatic scenarios supports the projection made by Rana et al. (2021) [3] in the Nepal Himalaya using ensemble modelling (eSDM) approaches. Shreshta et al. (2021) [53] proposed a contrasting viewpoint with HadGEM2-ES, CCSM4, and BCC-CSM1-1 simulation model in a different study. Of the three models, BCC-CSM1-1 and CCSM4 simulations anticipated the species to lose 61-85% and 71-98% of its niches by 2070, whereas, with HadGEM2-ES, the species will lose all preferred niches by 2070 in RCP4.5 and RCP8.5 scenarios. In such conditions, parallel studies on community structure, biotic interactions, plant-pollinator or plant-mycorrhizal fungi interactions, population status, and habitat characteristics are also necessary. Hence, an integrated model with both the elements (biotic and abiotic) and their cumulative effects needs further investigation. This is notable because they play an important role in the persistence of D. hatagirea populations [52].
Along with climate change, irrational harvesting practices of MAPs as a result of increased market demand for herbal medicine and lack of knowledge/ awareness on sustainable harvesting led to serious habitat degradations in the Himalayan region [49]. Similarly, changes in LULC, anthropogenic interference (i.e., habitat degradation and fragmentation, unscientific collection, pre-mature tuber overharvesting, grazing), interspecies competition (i.e., mainly habitat engulfed by Persicaria wallichii), and poor seed germination has seriously impacted the relative distribution of D. hatagirea in the study area [49]. The present study provides an overview of habitat suitability and likely changes with projected climate scenarios in conjunction with changes in the species geomorphologic profile. In this context, the information provided in the present study could be beneficial in various conservation initiatives at the local and regional levels in West Himalaya. The study conducted has some limitations as the model developed to forecast the fundamental niche of the species rather than the realized niche. The species realized niche might not be the same as predicted in our model prediction results. Another limitation is that this study modeled the habitat suitability of D. hatagirea at different time scales (2050 & 2070) based on abiotic factors but did not consider biotic factors such as plant-mycorrhizal association or plant-pollinator interaction. As a result, potentially suitable area of the species might be overestimated, as biotic interaction, especially mycorrhizal association, play a crucial role in plant establishment in the family Orchidaceae [53,73].

Conclusion
Dactylorhiza hatagirea, an endemic high-value Himalayan medicinal species, is under great stress and needs urgent attention. The perturbance of climate change has added extra pressure on the species allied with high degrees of anthropogenic stress. Therefore, the present study provided a clear overview of the habitat suitability for D. hatagirea and predicted the potential impacts of future climate on their distribution in Western Himalaya. The species potential distributions are explained mainly by elevation, precipitation of the coldest quarter, soil moisture, LULC, and mean annual temperature. The study reveals that D. hatagirea has approximately about 131 sq. km as its fundamental niche with high (>0.75) habitat suitability, which corroborates to about 0.24% of the total geographic area of the state. Under future projections RCP 4.5 and RCP 8.5, species distribution is projected to be very similar to the current distribution, except shifting species in a northward direction to higher elevation is expected. Furthermore, most predicted potential habitats fall within areas with anthropogenic encroachment leading to habitat degradation, unscientific and destructive harvesting practices, grazing, etc. Considering these, it is anticipated that the species may lose much of the habitat due to anthropogenic activities, while climate change impact will invariably be more at comparatively lower altitudes. Thus, it is imperative to undertake appropriate conservation steps and interventions like reintroduction/ augmentation programs. Besides these, the dependent communities awareness and sensitization curriculum holds utmost importance towards sustainable utilization. The model generated suitability maps can be of significant help to various policy prescribing National agencies, i.e., National Medicinal Plant Board (NMPB), Ministry of Environment Forest and Climate Change (MOEF &CC), Department of Science & Technology (DST). At the State level, the State Government has control over Forest Department and NGOs, thus ensuring and safeguarding the overall health of our forests and meadows.
Supporting information S1 Table. Multi-collinearity test by using cross-correlations (Pearson correlation coefficients, r) among environmental variables using ENM Tools.